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A Convex Feasibility Approach to Anytime Model Predictive Control 

Alberto Bemporad, Daniele Bemardini, Panagiotis Patrinos 


Abstract —This paper proposes to decouple performance opti¬ 
mization and enforcement of asymptotic convergence in Model 
Predictive Control (MPC) so that convergence to a given terminal 
set is achieved independently of how much performance is opti¬ 
mized at each sampling step. By embedding an explicit decreasing 
condition in the MPC constraints and thanks to a novel and very 
easy-to-implement convex feasibility solver proposed in the paper, 
it is possible to run an outer performance optimization algorithm 
on top of the feasibility solver and optimize for an amount of time 
that depends on the available CPU resources within the current 
sampling step (possibly going open-loop at a given sampling step 
in the extreme case no resources are available) and still guarantee 
convergence to the terminal set. While the MPC setup and the 
solver proposed in the paper can deal with quite general classes of 
functions, we highlight the synthesis method and show numerical 
results in case of linear MPC and ellipsoidal and polyhedral 
terminal sets. 


1. Introduction 

Model Predictive Control (MPC) is a well known advanced 
control approach in industry for its capability of optimizing 
closed-loop performance subject to operating constraints on 
input and output variables 0-0. In recent years, MPC has 
become very attractive also in fast-sampling applications with 
stringent real-time requirements, such as those arising in the 
automotive and aerospace industries. Such requirements posed 
a research challenge for developing optimization algorithms, 
and in particular Quadratic Programming (QP) solvers, that 
enable the use of MPC in commercial products. In particular, 
an embedded optimization solver must be fast, simple to code 
and test, require little memory, and have good worst-case 
estimates of its execution time. 

To cope with such requirements, multiparametric QP was 
proposed in 0 to pre-solve the QP off-line, therefore con¬ 
verting the MPC law into a continuous and piecewise affine 
function of the state vector. The main drawback of explicit 
MPC is that it is limited to relatively small problems and to 
linear time-invariant (LTI) systems. 

On-line optimization methods like active-set methods Q- 
Q, interior-point methods |[8|-|T0t, and dual piecewise 
smooth Newton methods m can be very effective in speed, 
but their worst-case CPU time can be hard to estimate in a 
non-conservative way. For accelerated dual gradient-projection 
methods p^ , good bounds on the worst-case execution time 
were provided (ID, (E), although the methods act on the dual 
QP problem, and therefore can lead to infeasible solutions if 
the execution is interrupted. 

On the other hand, in real-time control platforms the time 
allotted for the MPC controller to run is often not enough to 
cover the worst-case execution time, and other higher-priority 
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tasks may even preempt its full execution. Driven by such real¬ 
time constraints, anytime control algorithms were developed 
in GD with the idea of storing a set of control laws, each 
one of different complexity and closed-loop performance, and 
execute the one whose complexity is compatible with the 
current available CPU resources. 


In this paper we propose instead an MPC approach based 
on anytime optimization, with a novel convex optimization 
algorithm that recursively finds feasible solutions of decreasing 
level of suboptimality, depending on the computation power 
available within the sampling step. We first prove a rather gen¬ 
eral recursive feasibility and convergence result of MPC based 
on stability constraints that artificially impose a certain Lya¬ 
punov function to be decreasing |T^ , fT7| , where this function 
might be totally decoupled from the value function typically 
considered for assessing asymptotic convergence (H, GD- 
Moreover, in this paper we focus on convergence to a set 
around an equilibrium rather than to an equilibrium state, 
as in practical applications is often sufficient to track a set- 
point within a given tolerance |T9| . In addition, contrarily 
to |T0| that guarantees feasibility in real-time through a warm¬ 
starting technique in combination with robust MPC design, 
we provide an approach based on an original method for 
solving unconstrained problems, which is used for finding a 
feasible solution to a set of convex constraints. This method 
is very efficient in speed and easy to code, and can be run 
multiple times to approach an optimal solution, depending on 
the available CPU time. 


The paper is organized as follows. Section defines the 
main MPC setup and states recursive feasibility and conver¬ 
gence results. Section |nl| presents the new convex feasibility 
and optimization algorithm setup and shows its properties. 
Section |IV] proposes two ways of synthesizing a proper ter¬ 
minal set and stability constraints for linear systems subject 
to linear constraints on inputs and outputs, and Section |V] 
shows numerical evidence of the advantages of the proposed 
approach. 


A. Notation 

The sets of real and nonnegative integer numbers are 
denoted by M, N, respectively. For a vector x G denotes 

the i-th entry of x, and the expression x > 0 means that 
> 0, Vi = 1,..., n. For a matrix A G Ai denotes 

the i-th row of A, and A > 0 positive definiteness of A. 
Given a scalar r, denotes max{r, 0}; for a vector x G R"^, 
is the vector whose coordinates are = max{xi,0}, 

Vi = 1 ..., m. 



11. Feasibility-based MFC 
Consider the problem of steering the system 

x(t + 1) = a{x{t),u{t)) (1) 

to a target set 5 C while satisfying the constraints 

g{x{t),u{t)) < 0 ( 2 ) 

for all t G N, where x G u G a : ^ M and 

g : ^ M. We represent the target set as 

S = {xeR^ : f{x) < 0} (3) 

where f : ^ R and S is constrained controlled invariant 

with respect to in accordance with the following definition. 

Definition 7.' A set 5 C is constrained controlled 
invariant if for all x G 5 there exists u G R^ such that 
g{x,u) < 0, a{x,u) G S. 

Note that in (|^-([^ we are assuming scalar functions /, g 
without loss of generality. In fact, for any vector function 
h : R'^^ component-wise constraint h{x) < 0 is 

equivalent to the constraint h{x) = max^=i^...^^2 hi{x) < 0. 
Note also that, given any set 5 C a corresponding function 
/ : ^ M satisfying © can be defined as f{x) = 75 (x) — 1 , 

where 75 (x) = inf{/i : g > 0, x G gS} is the Minkowski 
function of S. 

To solve the stated control problem, we consider the fol¬ 


lowing MFC formulation 

N-l 

min £n{xn)+ X] 4 {xk, Uk ) (4a) 

fe =0 

s.t. xq = x{t) (4b) 

^/c+i = a{xk,Uk), /c = 0,1,. .., - 1 (4c) 

g{xk,Uk) <0, A: = 0, 1,..., - 1 (4d) 

/(xat) < 0 (4e) 

N-l 

^ f{xk)+ < (t){t - 1) (4f) 

k=l 


where 4 : R^+^ ^ R, k = 0,..., N - 1, In : ^ M 

are stage and terminal costs, respectively, and 0(7) is a given 
scalar, chosen in accordance with the following theorem. 

Theorem 2: Let u{t) = Uq be the control input applied to 
the process ([T]), where is any feasible solution of 

problem <0 at time f, and the quantity 

AT-l 

- 1 ) = /( 4 “')+ ( 5 ) 

k=l 

is constructed from the previous feasible solution , x]^^ 
of problem 0 at time f — 1, for alH G N. If the set S defined 
in ^ is constrained controlled invariant and problem Q is 
feasible at time 7 = 0 for the initial state x(0) and some value 
0(—1), then it is feasible at all time 7 G N and x{t) ^ S for 
t ^ oo. 

Proof: Let be any feasible solution of prob¬ 
lem (0 chosen at time 7, and let be the corresponding 

state trajectory. Consider the following candidate feasible 
solution {uk}^SQ at time 7 -f 1, where uq = u\, ui = 


..., un -2 = h]sf-i such that g{x)^,U]sf-i) < 

0, which exists by constrained controlled invariance of S. 
Let the state trajectory corresponding to {uk}^SQ, 

with Xq = x{t + 1). By construction, Xk = for 

all /c = 0, ...,A7 — 1, and hence g{xk^Uk) < 0 for all 
k = 0,..., A7 — 2, f{xN-i) < 0. Moreover, by the choice 
of ujsf-i we have g{x]sf-i,U]sf-i) < 0 and f{x]sf) < 0. Since 
u{t) = uo{t), and hence x\ = x(7 + 1), we have 

N-l N-l 

fi^k)+ = /(4)+ = </>(*) - fi^it+ 1))+ (6) 

k=l k=2 

SO that, since /(x(7+l))+ > 0, also the stability constraint 
is satisfied. Therefore, problem 0 admits a feasible solution 
at time 7 + 1, and because of 0. whatever is the choice of 
at time 7 + 1, we have 

<p{t + l) <4){t) - f{x{t + l))+ (7) 

This proves that limt^oo 0(7) exists, as 0 is a monotonically 
decreasing sequence and lower-bounded by zero, which in 
turns implies by 0 that f{x{t))+ = 0. If by 

contradiction we assume that x{t) -ft S for 7 ^ cx), then 
a subsequence 7/^ G N, G N, and a scalar (5 > 0 exist such 
that f{x{th)) > S, \/h G N, or equivalently f{x{th))+ > 
which contradicts limt^oo /(^(7))+ = 0. ■ 

Note that the convergence result of Theorem does not 
involve at all the cost function ( |4a| ). Of course, the transient 
behaviour of the system depends on how close to optimality 
are the chosen feasible solutions K}£o'of0. 

While the result of Theorem 0 does not make any as¬ 
sumption on the properties of functions a, /, g, Ik (except for 
constrained controlled invariance of <S), from now on we will 
restrict our attention to convex functions /, g, Ik and linear 
functions a, i.e., linear models 

x{t + 1 ) = Ax{t) + Bu{t) ( 8 ) 

in order to solve problem 0 effectively, using the novel 
algorithm proposed in the next section. 

III. Convex Feasibility Algorithm 
C onsider the following feasibility problem: 

Find X G C = {x G : fi{x) <0, 7 = 1,...,m}, (9) 

where : R'^ ^ M, 7 = 1,..., m are convex, twice continu¬ 
ously differentiable functions. Froblem 0 can be reformulated 
as the following unconstrained minimization problem: 

min F(x) = ff{x)+\\l (10) 

where f{x) = [fi{x) ... Nix)]. 

Proposition 3: If C 7 ^ 0, then inf F = 0 and arg min F = 

C. 

Clearly, if inf F > 0 then C is empty. 

Lemma 4: Function F{x) = ^||/(x )+||2 is convex and 
continuously differentiable with 

m 

wF{x) = vf{xyfix)+ = y2Mx)+^M^) (11) 


Proof: Function F can be written as 

^ m ^ m 

= -^(max{/i(a;),0})2 = , 

where R 3 z q{z) = z‘^ and R 3 z \-3 'ipi{x) = 
max{/^(x), 0}, q o is convex since q is convex and non¬ 
decreasing when its argument is nonnegative and fii is convex 
(as the pointwise maximum of the convex function fi and 
0) and nonnegative. Therefore F is convex as the sum of 
convex functions. On the other hand, F{x) = 
where R 3 z \-3 ip{z) = Since cp is continuously 

differentiable with continuous differentiability 

of F and formula (Tl\ readily follow. ■ 

Our ultimate goal is to devise Newton-like methods for solving 
the unconstrained problem ( p^ and thus ([^. Function F is 
but not however its gradient \/F{x) is a piecewise smooth 
mapping, in the sense that for any x eBF we have 

VF(x) G {\/Fj{x)}iex. (12) 

where I is the collection of all subsets I C {!,... ,m} for 
which there exists a x G such that fi{x) > 0, i G / and 
fi{x) < 0, i ^ /, whereas VF/(x) = fi{x)\/fi{x). The 

pieces of VF are smooth with Jacobian given by 

V^Fi{x) = Y,V fi{x)V fi{x)' + fi{x)V^fi{x). (13) 

iei 

For any x G let I{x) = {i G {1, • • • ,m} : fi{x) > 0}. 
Then the matrix will serve as a generalized Hessian 

of F at X, furnishing a second-order approximation similar to 
the one provided by the classical Hessian for functions. 
Another idea, stemming from Gauss-Newton methods for 
solving least-squares problems, would be to use as a gen¬ 
eralized Hessian the matrix fi{x)\/fi{xy, which 

results by omitting second-order terms fi{x)fi{x) from 
This choice saves us from computing the Hessians 
of fi, i ^ I{x) (notice that in case of fi{x) = a'-x — bi this 
makes no difference). However, this is a good choice only if 
we know that C is nonempty. In this case for any x G C 
we have fi{x) = 0 for i G I{x) and the term fi{x)fi{x) 
vanishes. 

Algorithm is a regularized piecewise smooth Newton 
method with line search. Its convergence properties can be 
inferred as a special case of | [20| , Specifically, every accu¬ 
mulation point of the sequence generated is a stationary point 
of F, and if is nonsingular then the convergence 

rate is quadratic. 

A. Feasibility-based optimization 

Problems involving constrained minimization of a con¬ 

vex function can be attacked by solving a sequence of feasi¬ 
bility problems. Specifically consider the problem 

min /o(x) (16a) 

s.t. X e C (16b) 

where C C is a closed convex set described as in (|^. Let 
/* = infxec fo{x) and X* = argmin^^ec fo{x). We assume 


Algorithm 1 : [empty, m]=PSN_FEAS(C) 

Input: a G (0,1/2), ( G (0,1), x^ G k = 0 

1 if F{x^) = 0 then 

2 I empty^false, x 3- x^\ exit 

3 else if ||VF(m^)|| = 0 then 

4 I empty=true; exit 

5 end 

6 Compute G that solves 

{V^Fik (x^) + 6^I)d = - VF(m^), (14) 

where = {i G [m] \ fi{x^) > 0}, 6^ = C||VF(m^)||. 

7 Compute = max{2“* | i = 0,1, 2,...} such that 

F{x^ + r^d^) < F{x^) + aTk\/F{x^yd\ (15) 

8 ^ x^ F T^df 

9 k 3- k and go to Step 1. 


that the set of optimal solutions of is nonempty. We have 
that 

/* = inf{t : x e S{t)}, X* = S{f*) 

where 

S{t) = {x€C: foix)<t} (17) 

is the lower level set of /o over C. Obviously we have that 
f > /★ if and only if S{t) is nonempty. This suggests that 
we can test whether a given t is smaller or larger than by 
solving the feasibility problem of finding x e S{t) 

m 

Fit) = 7^3;) = l{fQ{x)-t)l + iy2(fii^)F (18) 

i=l 

using Algorithm The following proposition, whose proof 
is omitted here for lack of space, proves some interesting 
properties enjoyed by function p. 

Proposition 5: Assume problem admits an optimal 
solution and let function p be defined as in ( p^ . 

(i) p is real-valued with p(t) > 0 for f < /*, whereas 
pit) = 0 for t > /*, 

(ii) p is convex and continuously differentiable with 

= Vt7(i, Xt) = -{fo{xt) - t)+, (19) 

where xt G argmiiijjeEn 7 (x,t). 

(iii) < 0 for t < /*, whereas = 0 for t > /*. 

Proposition |5]|^ shows that /* is the left endpoint of the 
halfiine {tGM: pit) = 0}. Therefore problem ( p^ has been 
reduced to finding the leftmost zero of the one-dimensional, 
monotone decreasing function : M ^ M. One way to find 
/* is to apply bisection to p. Starting from an initial closed 
interval with t- < we pick the midpoint 

t = (t_ + t+)/2 and try to determine if the level set Sit) 
is nonempty, by solving © using Algorithm (of course 
we could apply any other algorithm for unconstrained 
optimization). If pit) = 0 then Sit) is nonempty and this 
means that the corresponding Xt G argmin^^^Mn 7 (t, x) is 









feasible for (GUf, therefore t > and the new interval is 

reduced to In the case where ip{t) > 0, S{t) is empty, 

meaning t < so the new interval becomes 

1) Strengthening the lower bound: Overall, the bisection 
algorithm maintains a lower and upper bound for /*. Since 
the interval is halved at every bisection step, we obtain the 
standard linear convergence for bisection, that is, the algorithm 


stops after at most 


log2 


steps, where e > 0 is the 


desired optimality threshold. However, with almost no extra 
effort we can do much better in practice. 

Suppose that t < /*. The optimality condition for prob¬ 
lem is = 0 or 


{h{xt) - t)+Vfo{xt) + ^(/i(rct))+V/i(xt) = 0. (20) 

i=l 

From Proposition [5]|iii|), we have that f{xt) > t. Dividing by 
fo{xt) -t in (|2^ and letting m i = 1, • • •, m, 

we obtain 

m 

Vfo{xt) + Y,f^i^M^t)=0. 

i=l 


Since /i > 0, it follows that /i is a dual feasible vector, 
therefore 

m 

foixt) + < /*• (21) 

i=l 


Hence tn = fo{xt) + provides a lower bound 

on /^. Since to - t = fo{xt) + jZ'iLi > 0, to is 

indeed a tighter lower bound to /* than t. 

2) Strengthening the upper bound: When t > , due 

to [Sjl^, we have that f{xt) < t, Xt ^ arg mina^^i^r*. 7 (x, t). 
Therefore, if ip{t) = 0, Xt is a feasible vector for problem ( p^ 
and f{xt) is a tighter upper bound to /* than t. 

In fact, we can do even better. At every step of bisection we 
have at our disposal two vectors xp and xj, corresponding to 
the upper and lower bounds on /*, respectively. Vector xp is 
feasible, i.e., fi{xp) < 0,i = 1,..., m, while xj is infeasible, 
i.e., fi{xi) > 0 for at least one i and fo{xi) < fo{xp) (this 
follows directly from ([21])). Invoking a result by Bertsekas p2{ 
Proposition 2] we have that 

/* < |^^/o(a:F) + ^ Mxf), 

where T = inf {7 > 0 | fi{xi) < -'jfiixp), i = 
1,... ,m} Proposition 2]. The bound is nontrivial, i.e., 
the rightmost inequality is strict, when T < oc, which holds 
if and only if fi{xi) < 0 for all i with fi{xp) = 0. In that 
case we have T = inax|i|j.(^^)<o} ' 

The proposed improved bisection method is summarized in 
Algorithm 

3) Equality constraints: The approach can be immediately 
extended to handle linear equality constraints 


c'^x = dj, j = l,...,me 

in Q, by simply adding the term | P,^li{CjX — dj)^ in 
and to ip{t) in ( p^ , respectively. 


Algorithm 2: PSN OPT(/o, C) 


Input: accuracy e, xp e C, = fo{xp), xj ^ C, 
t- = Mxi) with 

Output: Xp e C with f{xp) — f'^ < e 

while —1_ > e do 

t i — (t_ + tp)l2 

Call [empty,Xt]=PSN_FEAS(5'(t)) (Algorithm 
if empty=true then 

XI ^ Xt 

t- ^ fo{xi) -I- l^ifi{xi), where 

_ {fi{xi)) + 
fo(xi)-t 

if fi{xi) < 0 for all i with fi{xp) = 0 then 
*+ ^ i^Mxf) + fyfoixi), where 

end 


1 
2 

3 

4 

5 

6 

7 

8 

9 
10 
11 
12 

13 

14 

15 

16 

17 

18 end 


else 


Xp ^ X 

if fi{xi) < 0 for all i with fi{xp) = 0 then 
t+ -s- ryfoixp) + r^fo{xi), where 


else 


■p _ vvno'iF fijxi) 

r — maxp|j7a.^)<o} -f^xp) 


t+ ^ Mxf) 


end 


end 


4) Determining initial upper and lower bounds: To deter¬ 
mine an initial upper bound to /* we can solve the feasibil¬ 
ity problem (|^ using Algorithm Determining a lower bound 
is a more delicate issue. If fo{x) = {l/2)x'Qx -P q'x where 
Q is symmetric positive definite, we can simply determine a 
lower bound on /* by computing the unconstained minimum 
X = —Q~^q. In general, if /o is convex and coercive we can 
find a X G such that V/(x) = 0. The nonlinear system 
can be solved by Algorithm In the case of a quadratic 
program with the cost having a positive semidefinite Hessian a 
lower bound to f'^ can be determined by solving the following 
convex feasibility problem 

m 

find X G /i G s.t.\/fo{x) + ^ fi{x) = 0, /i > 0 

Then /i is a dual feasible solution and q = /o(^) + 
^ /* (even if strong duality does not hold). 
Another way to determine a lower bound for general convex 
problems is to find a dual feasible vector, corresponding to the 
primal feasible vector xp corresponding to 

m 

find n e M™, s.t.WfoixF) + 'E =0’ M > 0. 

Notice that the set of p satisfying the conditions above is 
polyhedral. 

5) Special cases: Algorithm and can be used to 
solve linear programs (LPs), quadratic programs (QPs), and 
quadratically-constrained quadratic programs (QCQPs). In this 




















case, the computations in (V2) and ( p^ require only matrix- 
vector products. 

B. Applicability to feasibility-based MFC 

Algorithm [^requires all the constraints in the inner feasibil¬ 
ity problem to be twice differentiable functions. In particular, 
constraint (|^ is not continuously differentiable because of the 
max operator, so the above algorithm cannot be directly ap¬ 
plied to solve 0- However, we can simply recast the problem 
by introducing N — 1 additional variables e/c,/c = 0,...,A^ — 1 
and replace with 

e* > f{xk), fc = 1, 2,..., - 1 (22a) 

e* > 0 (22b) 

N-1 

- f{x{t))+ (22c) 

k=l 

without altering feasibility and optimality of the solutions. 

Moreover, in case / is given as the max of convex functions 
fi : R, i = 1,..., n/, constraints ( |4^ and ( |22a| ) can 

be replaced by 

Mxn) <0, i (23a) 

^k>fi{xk), 2 = k = l,...,N -I (23b) 

Similarly, if g is given as the max of convex functions pi : 
^ M, i = 1,..., n^, ( |4d| ) can be replaced by 

9i{xk,Uk) <0, i = 1,..., /c = 0,..., - 1 (24) 

The case ^ik can also be dealt with by 

introducing "f2k=o additional variables aik, k = 0,..., N, 
and replacing ( |4a| ) with 

N riik 

min EE aik (25a) 

k=0i=l 

s.t. aik > ^ik(xk,Uk), i = 1,..., n^k, k = 0,... ,N - 1 

(25b) 

o'iN > ^iN{xk), i = 1,..., n£N (25c) 

In conclusion. Algorithm can be applied to solve 0 for any 
twice differentiable convex function fi, pi, iik. 

IV. Constrained tracking to a set 
C onsider an output vector 

y[t) = Cx{t) (26) 

associated with process 0, with ^ G M^, and a corresponding 
output reference r elZ CW, where 7^ is a polytope. Assume 
that the following linear system 

0 = AXr + BUr — Xr 

r = CXr 

admits a unique solution (xr^Ur) of steady-state state and 
input vectors for all r £ IZ, and assume that IZ is such that 

g{Xr^ Ur) < 0. 

We consider the problem of controlling ([8]) to a desired set 
At = {x e : S{x — Xr) < s}, around the equilibrium 


state Xr, where 5 > 0, s G , while satisfying the input 
constraints 

^min E u(t) ^ '^max (27) 

and the output constraints 

^min E yif) — ^max (28) 

with Umin < 0 < ^max, Vmin < 0 < ^^ax- In this CaSC, WC 
define the function g in ( |4d| ) as the convex and piecewise affine 
function 


g{x,u) = max {Gi {Ax-h Bu)G'^u — g^} (29) 


where G^ = 


0 

0 

c 

-c 


, G^ = 


I 

-I 

0 

0 


9^ = 


i^max 

?/max 
— ?/min 


, and q = 


2m-\-2p. An example of desired set AV is given by S = [ _^ ], 
s = [-e^rn]’ convergence to AV implies satisfying 

the constraint on the tracking error Cmin ^ 9 — r < Cmax 
asymptotically. 


A. Quadratic functions 

We consider the ellipsoidal terminal set S defined by 

f{x) = {x — Xr)' P{x — Xr) — Pr (30) 


where P = P' > 0 and pr are determined in accordance with 
the following theorem. 

Theorem 6: Let Q = Q' > 0, Q G Y G 

X = A' > 0, A G be the solution of the semidefinite 

program 


max (det Q) ^ 

Q 


s.t. 


> 0 


> 0 


{AQ + BY)' 

AQ + BY XQ 

X Y 
Y' Q 

i = 1, . . . ,TO 

Q {AQ + BYyC'i 

Ci{AQ + BY) 

Q Q'S'i 


v} 


SiQ sf 


i — 1,... ,nT 


where 0 < A < 1 is a given contractive factor, and 


(31a) 

(31b) 

(31c) 
(3 Id) 
(31e) 

(31f) 


min {riniax,i (32a) 

j = l,...,raR 

Vi ~ . min {t/maXji + {rj)i} (32b) 

j = l,...,nR 

for alH = 1,... , g, where are the vertices of IZ. Then, 

by setting K = YQ~^, P = Q~^, f{x) = {x — XrfPix — 
Xr) — pr, the set S = {x \ {x — Xr)'P{x — Xr) < pr} is 
constrained controlled invariant under the control law u = 
K{x — Xr) Y Ur in that x G 5 implies 


Ax Y Bu G S 

(33a) 

Win E K{x Xr) “h Ur E '^max 

(33b) 

^min E G{Ax Y Bu) Y ^max 

(33c) 

Sx < S 

(33d) 















for all r ^IZ, where 


in I - ^ - 1 > 1 (34a) 

,q+nT \ AiQA^ j 


K 


Umax Ur 

-K 


'Umin P Ur 

C{A + BK) 


^max T 

-C{A + BK) 


-ymin P r 

S 


5 


Proof: Let Ax = x — Xr, Au = u — Ur, Ay = y — r. 
Clearly Ax{t + 1) = Ax + Bu — Xr = AAx{t) + BAu{f) 
and Ay{t) = CAx{t), along with the constraints r^min —< 
Au(t^ ^ '^max ^min T "A CAx{f) ^ ^max The 
robust satisfaction of properties ( [33] ) with respect to r G 7^ for 
all X such that Ax'PAx < 1, under the control law Au{t) = 
KAx{t), follows by standard arguments from the inequality 
constraints in ( [ST] ) (see, e.g., |[23|). For a given r e IZ, the 
scalar defined by ( [^ provides the largest ellipsoid centered 
in Xr and defined by P such that ( [33] ) are satisfied, where 
Pr > 1. ■ 


B. Polyhedral terminal set S 

For a given under asymptotically stabilizing feedback con¬ 
trol law u(t) = K{x{t) — Xr) Ur, considcr the polyhedral 
terminal set S defined by 

f{x) = max{Ff-(x — Xr) — Ki} (35) 

where 5 = {x : P[{x — Xr) < Ff} = {Ax : HA < 
Ff} is a maximum admissible polyhedral invariant set p4| 
for the closed-loop system Ax{t -f 1) = {A BK)Ax{t) 
and with respect to the constraints A Ax < b^in, where A is 
defined in ( [34b| ), and 6niin is defined as in ( [34b| ) by replacing 
("^max Ur)i with rniUjf^]^^ {(U/j^iax Urj)A, ( 

with minj=i^. .^^^ {(- 'Umin + Ur.)i}, i = {y max 

r)i with {(2/max - rj)i} and (-//min + r)i with 

minj=i,...,„jj{(-2/mm + rj)i}, i = Clearly, the size 

of S depends on the size of At and on how large are the 
boxes {y e W : ^min < y < ^max} with respect to 7Z 
and {u G ^ u < r^max} with respect to the set 

{u G : u = Ur, r G P}. 


V. Simulation Results 


Consider the linear system described by the transfer function 


G{s) 


26(5 + 1) 
52 + 25 + 26 


(36) 


Model ( [ 3 ^ is converted to discrete-time by exact sampling 
plus zero-order hold with sampling time = 0.2 s, resulting 
in the state-space model with matrices A = [ 0.4424 ] ’ 

^ = [ 2.0623 ]’ ^ 1.9407], D = 0. The system is 

subject to the constraints 


(37) 

leading to defining g{x,u) as in ( [29] ). We setup the MFC 
problem <0 with A = 6, 

ik{x, u) = {x — XrYQix — Xr) p {u — Ur)'R{x — Xr) (38) 


Q = lOC'C, i? = 1, for all A: = 0,..., A - 1, and 

^n{x) = (x — Xr)' P{x — Xr) (39) 


where P is the solution of the discrete algebraic Riccati 
equation associated with A, B, Q, and R, and Xr = [ lAQsg ] r, 
Ur = r. The possible reference signals are restricted in 
the interval P = [—0.9,0.9], while the desired target set 
At = {x G : ||x — x^Hoo < 0.1}. We start from the initial 
condition x(0) = [q] and command the set-point r = 0.5. We 
consider the cumulated cost 

10 N-l 

J-Yl +(+) + ^fe(4>4)) 

t=0 k=0 


as a measure of closed-loop performance, where is 

the solution of problem ([^ chosen at time f, and {x^}^q the 
corresponding state trajectory. We consider two settings for 
function /: 

Case (i): the convex quadratic function as in ( [^ , where 
P = [^ 23 . 97^3 io¥ 4493 ] obtained by ([^ with A = 1, along 
with the terminal gain K = [ 0.1968 -0.2898], Tahle|i]shows the 
cumulated cost J for different maximum values of permitted 
CPU time to solve problem ([^. The corresponding trajectories 
are depicted in Figures [^ [2 


TABLE I 

Performance vs. allocated CPU time (ellipsoidal constraints) 


max CPU time (ms) 

Cumulated cost J 

unbounded 

8.7865 

60 

22.6572 

40 

26.4663 

20 

31.4528 


Case {ii)\ the convex piecewise affine function as in ( [35] ), 
where H = [l\ V o' '‘"f" = 0.1 [i i i i i i]' 

is the maximum A-contractive invariant set for the closed-loop 
system Ax(t + 1) = {Ap BKi^Q^)Ax{t), Ki^qr is the LQR 
gain associated with A, B, Q, and R, and A = 1, computed 
as described in [ [25| . Table [II| shows the cumulated cost J 
for different maximum values of permitted CPU time to solve 
problem 0 . The corresponding trajectories are depicted in 
Figures [^ [^ 

TABLE II 

Performance vs. allocated CPU time (polyhedral constraints) 


max CPU time (ms) 

Cumulated cost J 

unbounded 

9.0075 

10 

23.4491 

5 

26.4646 

1 

26.7551 


Finally, we compare the performance of the new solver 
described in Section [I^ (implemented in interpreted MATLAB 
code) against the commercial solver Gurobi 5.6.2 in 
solving QCQP and QP problems to optimality, and also 
to qpOASES for QP’s. We consider problems deriving 
from the MFC setup with ellipsoidal (QCQP) and polyhedral 
(QP) constraints described above, for an increasing prediction 
horizon N. The results are depicted in Figure [^ (QCQP case) 
and Figure [^ (QP case), respectively. 
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Fig. 1. System trajectories for varying values of available CPU time 
(solid=optimal), ellipsoidal constraints 


VI. Conclusions 

The contribution of this paper is twofold. From an optimiza¬ 
tion viewpoint, we have introduced a very efficient numerical 
solver that can solve convex feasibility and optimization 
problems, and that is at least an order of magnitude faster than 
commercial state-of-the-art (interior-point) solvers as the di¬ 
mension of the problem increases. By taking advantage of the 
way the solver computes an optimal solution via a sequence 
of convex feasibility problems, from a control viewpoint we 
proposed an MFC strategy for convergence to a terminal set 
that allows an anytime optimization philosophy, that is of 
improving the optimality of the control move with respect 
to a given performance specification only if CPU resources 
are available during the sampling interval. We believe that the 
approach has potential applications in embedded MPC systems 
where a large-enough time-slot for computations cannot be 
guaranteed a priori, a rather typical situation in multitask real¬ 
time systems. 
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Fig. 6. QP: Algorithmic (blue), Gurobi (red), qpOASES (green) 

























